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An accurate model of a vertical pillar quantum dot is described. The full three dimensional 
structure of the device containing the dot is taken into account and this leads to an effective two di- 
mensional model in which electrons move in the two lateral dimensions, the confinement is parabolic 
and the interaction potential is very different from the bare Coulomb potential. The potentials are 
found from the device structure and a few adjustable parameters. Numerically stable calculation 
procedures for the interaction potential are detailed and procedures for deriving parameter values 
from experimental addition energy and chemical potential data are described. The model is able to 
explain magnetic field dependent addition energy and chemical potential data for an individual dot 
to an accuracy of about 5%, the accuracy level needed to determine ground state quantum numbers 
from experimental transport data. Applications to excited state transport data are also described. 

PACS numbers: 73.63.Kv, 73.23.Hk 



I. INTRODUCTION 

Vertical pillar quantum dots exhibit an extremely wide 
range of interesting physics, particularly in the presence 
of a magnetic field perpendicular to the plane of the dot 
P, 13, Q • The Coulomb blockade allows single electrons 
to be added to the dot so the number of electrons can be 
set precisely, starting from one. For each number of elec- 
trons the ground state evolves through a series of tran- 
sitions as the magnetic field is increased. At zero mag- 
netic field the ground states are well described by Hund's 
first rule and maximum density droplet (MDD) states oc- 
cur when the field is increased to a few Tesla. Finally, 
fractional quantum Hall droplet and electron molecule 
states occur at high fields beyond about lOT. So, roughly 
speaking, there are four regimes but the detailed pic- 
ture is much richer. Many transitions occur in between 
each of the main regimes and each transition is accom- 
panied by abrupt changes in the spin and orbital angular 
momentum of the ground state. There is experimental 
P, i, i, i, H and theoretical 1, i, 0, i, Q evidence that 
these effects occur. However it is difficult to apply experi- 
mental techniques, like scanning tunnelling spectroscopy, 
to probe the corresponding quantum states directly. In- 
stead the ground state quantum numbers are found by 
comparing data from transport spectroscopy with calcu- 
lated results. This requires an accurate dot model that 
can be used to analyse data for an individual quantum 
dot. Development of an appropriate model is the purpose 
of the present work. 

Although quantum dots are often described as artifi- 
cial atoms, there is a very important difference between 
an artificial atom and a natural one: all natural atoms of 



the same isotope are identical but all quantum dots are 
different. This makes it quite challenging to model an 
individual dot accurately, even when its basic design pa- 
rameters are known, because manufacturing tolerances 
introduce fluctuations in the parameters and this can 
have a significant influence on dot states. The fact that 
the dot often consists of a small region embedded in a 
much larger device structure just adds to the difficulty of 
constructing an accurate model. The ideal model of an 
electrostatic quantum dot is a system in which electrons 
are constrained to move in the two dimensions parallel 
to the dot plane, are confined by a parabolic potential 
and interact via the Coulomb interaction. The proper- 
ties of this model have been examined extensively, see Q 
for a review. It gives a good qualitative description of 
dot behaviour but a more accurate model is needed for 
data analysis. One possibility is a device model in which 
the dot confining potential and in some cases the inter- 
action potential are computed, together with the quan- 
tum states, from a combined solution of the Poisson and 
Schrodinger equations. This approach has been used by 
many authors, O, [HI [13, [II] for example, and has given 
a great deal of insight into the generic behaviour of de- 
vices containing quantum dots. However, it is extremely 
difficult to use generic device models to deduce ground 
state quantum numbers from experimental data for an 
individual quantum dot. 

One of the difficulties encountered in analysing individ- 
ual dot data is that energies have to be computed to very 
high precision. The quantity measured experimentally is 
the gate voltage at which transport occurs. This is pro- 
portional to the chemical potential, /ijv which is the dif- 
ference of two dot energies: = En — En^i, where E^ 
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is the energy of the an A''-electron dot state. Further, the 
transport data is often presented as a difference between 
the gate voltages for two successive current peaks. This 
gives the addition energy, Ean = Ejq+i — 2Ejs[ + i?Af-i, 
the second difference of with respect to N . Be- 
cause energy differences are needed to compare theoret- 
ical results with experimental data the energies them- 
selves must be computed very accurately. Typically, the 
addition energy is around 3 meV while typical dot en- 
ergies are one or two orders of magnitude larger so that 
high precision values of the dot energy are needed for 
data analysis. However the potentials found in generic 
dot models depend on many parameters, dot dimensions, 
dopant densities etc, some of which are not known accu- 
rately. In principle, these parameters could be adjusted 
to fit the data but it is desirable to avoid fitting a large 
number of parameters. 



The present model lies between the parabolic confine- 
ment. Coulomb interaction model and the generic device 
models. The main idea is to develop a model that con- 
tains all the essential physics but depends on a small 
number of adjustable parameters. Three are used in 
the present work but it turns out that one of them is 
not very significant. The model is a parabolic confine- 
ment model, believed to be accurate for small numbers 
of electrons [ll|. Two of the parameters determine the 
parabolic potential but only one of them is significant. 
The interaction potential is determined from a simpli- 
fied device structure. The effects of finite thickness and 
screening on the interaction in this structure are included 
fully and a third model parameter is used to determine 
the screening. Essentially, the resulting model reduces to 
a two-dimensional, parabolic confinement model with an 
interaction potential that is very different from the bare 
Coulomb potential. The model gives an excellent de- 
scription of addition energy data, accurate to about 0.15 
meV. It has already been used to determine the ground 
state spin and orbital angular momentum of a pillar dot 
in the strong magnetic field regime [5|]. However ref. Q 
does not contain the detailed explanation of the model 
and analysis procedures that is given here. 



The model is described in section |TT] and the method 
used to calculate the quantum states is detailed in section 
mil Following this, parameter fitting procedures for ex- 
perimental addition energy and chemical potential data 
are discussed in section IIVI Then the application of the 
model to data in the low and high magnetic field regimes 
is described in section |V] and conclusions are stated in 
section rvTl Finally, two appendices deal with the calcula- 
tion of Green's functions and interaction matrix elements 
needed to find the quantum states. 




FIG. 1: Typical vertical pillar dot (top) and typical transport 
data (bottom). Labels indicate the Landau level filling factor 

V. 



II. DOT MODEL 

A. Vertical pillar dots 

A vertical pillar dot (Fig. [1] top frame) consists of a 
narrow pillar, typically 500 nm in diameter Q. The pillar 
contains a double barrier heterostructure (DBH) which 
provides confinement in the vertical (z) direction and is 
surrounded by a cylindrical metal gate which carries a 
negative charge that provides lateral {x,y) confinement. 
The pillar stands on a substrate which is heavily doped 
and there is also a heavily doped region at the top of 
the pillar. These heavily doped regions provide source 
and drain contacts that enable the transport properties 
of the device to be investigated. A unique advantage of 
the vertical geometry is that it is not sensitive to edge 
states which affect transport through lateral devices in 
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the high magnetic field regime. 

In a typical experiment, the source-drain current is 
measured as a function of the gate voltage, Vg, source- 
drain voltage, Vsd and magnetic field, B parallel to the 
pillar. Typical results are shown in the bottom frame 
of Fig. [TJ This image shows the condition for an elec- 
tron to enter the dot at fixed source drain bias. It is 
well known that transport through the device is Coulomb 
blockaded except at certain values of Vg and Vsd [!]■ In 
Fig. [1] the current is colour coded so that dark regions 
correspond to the largest current. The A'^th dark curve 
shows how the gate voltage needed for the iVth electron 
to enter the dot depends on the magnetic field. The low- 
est curve in the figure corresponds to iV = 1. 

The exact condition for electron transport is 



eV 



sd 



> 



^iN - ea{VgN)VgM > tJ.c 



eVs 



sd 



(1) 



where /Zc is the contact chemical potential, VgN is the 
gate voltage at which the A^th electron enters the dot 
and a is an electrostatic leverage factor. At zero source 
drain bias this reduces to 



ea{VgN)VgN = ^ 



(2) 



so that the gate voltage at which transport occurs is a 
measure of the dot chemical potential relative to the con- 
tact. For transport data in the form of a gate voltage 
difference, VgN — VgN-i is proportional to the addition 
energy, Ean = E^+i — '2.En + En^i if a and fic are 
independent of N. Hence the first and second differences 
of the dot energies are needed to compare theory with 
experiment. 



B. Overview of model 

The dot model takes account of vertical confinement, 
lateral confinement and interactions between dot elec- 
trons. The quantum well between the double barriers is 
relatively narrow (12 nm) so the electrons are confined 
to the ground state of the vertical motion. This leads to 
a quasi-two-dimensional model in which the interaction 
is modified by the vertical confinement, physically the 
effect is to smear out the Coulomb singularity. 

The lateral confinement is generated by the cylindrical 
side gate and in principle the lateral confinement poten- 
tial ma y b e found by solving the Poisson equation, refs. 
[13, [0,112,1131, fo'^ example. However in the present case, 
N is small so the spatial extent of the dot state (< 50 
nm) is very small compared to the pillar diameter (w 500 
nm). Therefore the only part of the confining potential 
that has a significant effect on the dot electrons is the 
part near the minimum which is parabolic. In addition, 
observation of shell structure in the device considered 
here [^ shows that the device has a high degree of circu- 
lar symmetry. The lateral confining potential is therefore 
taken to be Vc = Vo+m*ujQr'^ /2, where r is the lateral ra- 
dial co-ordinate and luq is found from huJo = a + b{N — 1) 



where a and b are fitting parameters. The constant Vq is 
not needed for the present analysis because the addition 
energy is independent of Vq and the chemical potential 
is also independent of Vq when /ijv taken relative to fii. 
Previous work has shown that the parabolic approxima- 
tion is consistent with experiment [1| and accurate for 
small numbers of electrons that is electrons in the 
first and second shells at B = 0. 

The screening is mainly caused by heavily doped con- 
tacts that are relatively close to the dot, around 10 nm 
from the well. This is close enough for the contact re- 
gions to have a significant screening effect on the inter- 
action between the dot electrons. This is treated in the 
Thomas-Fermi approximation and the screening length 
in the contact region is taken to be a fitting parameter. 
Full details of the present approach to screening and fi- 
nite thickness are given in the following two sections. 



C. Finite thickness 

The dot is three dimensional but when the electrons 
are confined to the vertical ground state their motion is 
quasi- two-dimensional. A variational method is used to 
find the resulting effective Hamiltonian. In principle this 
needs two steps. The envelope function for the verti- 
cal confinement should be found first and then the spin 
splitting should be found from this envelope function [lj| . 
However the second step is not needed here because the 
experimental value of the effective ^-factor for the present 
device is known. Except for the spin splitting terms, the 
3D effective mass Hamiltonian is 
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2m*[z) 2m*(z) 



TT^ + Vein) + VbOiz^) 

(3) 



where tt — p|| -f eA, p|| is the lateral momentum, A is the 
magnetic vector potential, pz is the vertical momentum, 
Vb is the barrier height, 0(z) is a step function, Q{z) = 
1 when \z\ > w/2, Q(z) = otherwise, w is the well 
width, m*(z) = when z is in the well and m*{z) = 
ml when z is in the barrier. U is the electron-electron 
interaction potential and r = ix,y). The variational trial 
function for the electron states is taken to have the form 

— Ilix{zi)^{ri, ...tn), where $ is antisymmetric and 
includes spin functions. 

Equations for <I> and x are found by minimising 
(^'|i/|'I') subject to the constraint that $ and x are nor- 
malised. This leads to 



i?lll$) = A|$), 



E\\{ml) - E^\im*J 
N 



e 



lx)=A'|x), 



(4) 



where A and A' are eigenvalues, mj^, and ml are effective 
masses in the well and barrier respectively and Hu = 
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PwH\i {m^) + (1 -p,„)iJ|| {ml). Here is the probabihty 

of finding and electron in the well, = /^^yj ^^(^)'^'^ 
and the parallel and perpendicular Hamiltonians are de- 
fined by 
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+ ^E/ x\z)x\z')U{v,~v,,z,z')dzdz', 



2m* (z) 



Pz + VbOiz). 



(5) 



The parallel energies in the second of Eqs. ^ are ex- 
pectation values of i/||: E\\{m) = (m)|$). After 
Eqs. ([U have been solved the total energy is found from 
the expectation value of the full Hamiltonian: 



E^N{x\H^\x) + {mn\<^). 



(6) 



In principle, the coupled eigenvalue problem defined 
by Eqs. Q can be solved iteratively to arbitrary accu- 
racy but in practice this is not necessary. The energy 
difference term [E\\{ml) — E\\{ml^)]/N in the second of 
Eqs. ([3]) in the present case is around 1% of the barrier 
height. Hence the effect of this term is small and it is 
sufficient to find x from the zeroth-order approximation 
in which the energy difference term is neglected and then 
use the zeroth-order x to find the averaged Hamiltonian 
H\\. The averaging has two physical effects. First, it 
changes the effective mass. The dot considered here is 
made from InGaAs containing 5% In, and the In reduces 
the mass. The barrier is made from AlGaAs containing 
22% Al and this increases the mass. The net effect is to 
reduce the mass to 0.0653mo, slightly less than the GaAs 
value. The second effect of the averaging is to smear out 
the Coulomb singularity in the effective two-dimensional 
interaction. This results from the integral over z in the 
first of Eqs. (O. The averaging procedure also applies to 
systems with spin-orbit coupling. In this case, it is found 
that the Dresselhaus spin-orbit coupling parameter has 
to be averaged in a way similar to the effective mass. 
This has been used for studies of spin relaxation in the 
present device (l5| . 

The final step in determining the effective Hamiltonian 
is to find the spin splitting. The effective g-factor, 
is determined from the experimentally observed Zeeman 
splitting d in the same device that is used for the data 
analysis performed here. \g*\ = 0.3 for the dot when 
B > 10 T, compared with \g*\ = 0.44 for bulk GaAs. 
The reduced value for the dot is consistent with the ef- 
fect of non-parabolicity The non-parabolicity also 
makes the effective ^-factor depend on magnetic field and 
Landau level index but the resulting corrections to 
/i M are probably smaller than experimental error. A con- 
stant g-factor is therefore used and the Zeeman energy 
is accounted for by adding g* ^bBSz to the total energy 
in Eq. [S] Here /is is the Bohr magneton and Sz is the z 



component of the total spin S. The effective ^-factor is 
assumed to be negative, as in bulk GaAs, so g* = —0.3. 



D. Screening 

The most important source of screening in the present 
device is the heavily doped contacts above and below the 
dot. The dielectric response of the disordered contact 
material is not well understood but can be treated ap- 
proximately in the static Thomas-Fermi approximation. 
This allows the screening length to be estimated from the 
density of states at the Fermi level, however this is not 
known accurately for the contact material. Hence the 
screening length is taken to be a fitting parameter which 
is determined from experimental data. 

The screened interaction between dot electrons is 
found from the electrostatic Green's function G(R, R') 
which satisfies 

- V-[(e(R)eoVG(R,R')] = -g^(R)e(R)eoG(R, R') 

+ ,5(R-R'), (7) 

where e(R) is the relative permittivity, go — 27r/As inside 
the contacts, go = outside the contacts and R = (r, z). 
Because a detailed theory of magnetic field dependent 
screening by the 3D contact material is not available, the 
screening length. As, is taken to be independent of the 
magnetic field. This approximation should be good when 
the screening length does not change significantly over 
the field range used in the experiments and it is assumed 
that this is the case. Eq. ([7]) is simplified by making use 
of the fact that the dot is an order of magnitude smaller 
than the pillar. This means that the Green's function 
is almost translationally invariant in the lateral direction 
and can be found to a good approximation by replacing 
the pillar by a stack of dielectric layers of infinite extent 
in the lateral directions, e and go then do not depend on 
r. So the Green's function is translationally invariant in 
the lateral directions and can be expressed as a Fourier 
transform: 

G(R, R') = ^ y G(g, z, z') exp(iq • (r - r'))dq. (8) 
The Fourier components of the Green's function satisfy 



d 
dz 



e{z)eQ — G{q,z,z') 



+ e(z)eo[g' + 9o(^)]G(g, z, z') = 5{z - z'). (9) 
The solution of this equation has the form 



G{q,z,z') = - 
G(q,z,z') 



W ' 
9iz)f{z') 
W 



z > z 



z < z\ 



(10) 



where / and g are solutions of the homogeneous equa- 
tion corresponding to Eq. ^ and W is their Wron- 
skian. These solutions are chosen to satisfy / — > when 
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TABLE I: Layer structure of the device modeL 



FIG. 2: Effective 2D interaction (solid line). Dashed lines 
show the limiting form and the 1/r Coulomb interaction. 



z — > +00 and 5 — > when z —00. Although the ana- 
lytic form of G is as given in Eq. (jlOp , severe numerical 
instabilities are encountered when this form is evaluated 
unless special precautions are taken. The problem is that 
the G has rapidly growing components which cause expo- 
nent overflow. Exactly the same mathematical problem 
is encountered in calculations of evanescent wave prop- 
agation in electron diffraction theory and the solution 
used here is a reflection matrix approach developed for 
calculations of reflection high energy electron diffraction 
(RHEED) This is detailed in appendix A. 

Once the Green's function has been found, the Fourier 
components of the effective two-dimensional interaction 
are obtained in the form of an integral: 



X^iz)x^iz')Giq,z,z')dzdz'. (11) 



The special form of the Green's function, Eq. (|10p . is 
used to simplify this integral so that it can be evaluated 
efficiently. Details are given appendix B. The form of the 
effective interaction in real space is very different from a 
pure Coulomb interaction. Because of the screening, the 
interaction is dipole-like at long range and decreases like 
1/r^ in the limit of large r. And the Coulomb singularity 
is removed by the finite thickness because the electrons 
are able to move out of the dot plane and avoid each 
other. The real-space effective interaction as a function 
of r is shown in Fig. [2l U2d{r) is calculated numerically 
from its Fourier transform U2d{q) ioi the layer structure 
shown in Fig. [T] (see also Table [l| and As = 10 nm, a 
typical screening length. The length unit is the 2D har- 
monic oscillator length parameter, A'^ — h/{2m*Vl) where 
= Wq -I- with ijjc the cyclotron frequency. The 
energy unit Ex — e^/(47reu,eoA) where e^, is the relative 
permittivity in the well that contains the dot. It is clear 
that U (r) is significantly different from the 1 /r Coulomb 
interaction and approaches the form in the limit of 
large r. 

The effect of the dielectric interfaces in the device 



Layer 



Composition 



Thickness 



Top contact 

Buffer 

Barrier 

Well 

Barrier 

Buffer 

Bottom contact 



n+ GaAs 
i GaAs 
Alo.22Gao.78As 
Ino.05Gao.95As 
Alo.22Gao.78As 
i GaAs 
n""" GaAs 



3 nm 
9 nm 
12 nm 
7.5 nm 
3 nm 



on the effective interaction is accounted for fully in the 
Green's function formalism. In addition, the dielectric 
interfaces have an effect on the vertical confinement po- 
tential. This results from the interaction of each electron 
with its own image charges and is also accounted for in 
the Green's function formalism. The image charge con- 
tribution to the vertical confinement is obtained from the 
non-singular part of the Green's function and in the 
present case the image charge contribution is 



Viiz) 



87re^eo Jo 



2qe^eQG{q, z,z) - — dqll2) 



where e = e^, when z is in the well, e — eb, the barrier 
relative permittivity, when z is in the barrier and e = 
(e^, -|- eb)/2 when z is at the interface between the well 
and the barrier. The image charge term leads to a small 
change in the vertical confinement which is taken into 
account via perturbation theory. 

Another source of screening in the present device is the 
metallic, cylindrical side gate. In principle, this breaks 
the translational invariance of the effective lateral inter- 
action however this effect is expected to be small. The 
magnitude of the effect can be estimated from the elec- 
trostatic Green's function for an infinite metallic cylinder 
[20I I . Numerical calculations of the Green's function sug- 
gest that the effect is < 1% within a region of radius 
^ 10 nm in the centre of the cylinder. The effect is 
clearly largest for electrons at the edge of the dot but 
the interaction between these electrons and electrons in 
the centre of the dot is suppressed by the screening ef- 
fect of the contacts. These considerations suggest that 
the effect of the gate is small for the small electron num- 
bers considered here {N < 5). However, the effect could 
be significant for larger electron droplets whose edge is 
closer to the metallic gate. 



E. Model Hamiltonian 

The effective Hamiltonian for the two-dimensional dot 
model is 



i 

+ lT.U2dir^-r,), (13) 
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TABLE II: Material parameters for the device model. 



Parameter Material Value 

m*/mo AUGai_^As 0.067 + 0.083x 

m*/mo InyGai_j;As 0.067 - 0.045y 

Vb (meV) AUGai-^As / InyGai-^As 823a:: + 640y 

£(, AUGai_^As 12.7-3.122; 

e„ IriyGai-yAs 12.7 + 2y 



where Ez — (xI^-l|x) + AV/ is the perpendicular en- 



ergy, AV/ is the perturbation caused by the image charge 
term. The averaged effective mass m* is found from 
1/m* = pw/mlj + {l—pw)/'ml and the effective factor 
is determined experimentally (Section III C[) . The layer 
structure of the device is detailed in Table [J and the ma- 
terial parameters used for the present calculations are 
detailed in Table HJ 



III. CALCULATION OF EIGENSTATES 

The effective Hamiltonian H^ff is very similar to the 
Hamiltonian for a two dimensional parabolic dot with 
a Coulomb interaction. Hence the eigenstates of i?e// 
are found in the usual way by numerical diagonalization 
in a Fock-Darwin basis IS"]. The only difference between 
the usual calculation and the present one is that matrix 
elements of the effective interaction U2d have to be com- 
puted. In the case of the Coulomb interaction these ma- 
trix elements are normally found from the Fourier trans- 
form of the interaction and the calculation involves an in- 
tegral over g Q . In the present case the Fourier transform 
approach is also used. U^diq) has the form F{q)V2d{<l) 
where V2d{<l) is the Fourier transform of the Coulomb 
interaction, e^/(47reu)eor), in two dimensions and F{q) 
is a form factor that is computed numerically from the 
Green's function. The only difference between the stan- 
dard treatment of the interaction in a 2D dot and the 
present treatment is the appearance of the form factor 
in the q integral. This integral is done numerically by 
Romberg integration. 

The Fock-Darwin states are labelled by an angular mo- 
mentum quantum number / and a radial quantum num- 
ber n. The many-electron basis used for the diagonal- 
ization consists of Slater determinants formed from these 
states. It is important to minimise the size of this basis as 
expensive calculations have to be performed repeatedly 
to fit the model parameters. The Hamiltonian is block 
diagonalized according to the value of the total orbital 
angular momentum J. For iV < 4 the basis for each J 
is formed from Fock-Darwin states with n < 3 and all I 
values compatible with the required J value. For N = 5 
the size of the basis is limited by making use of the fact 
that the radial excitation is an inter-Landau level exci- 
tation and hence has large energy in a strong magnetic 
field. This enables the maximum n quantum number to 
be reduced as a function of magnetic field. In the present 



case the maximum value of n is taken to be the integer 
part of 6.9 — 0.28i3 which corresponds to a maximum n 
of 6 at B = T and 2 at S = 14 T. All Slater determi- 
nants compatible with this n value are then constructed 
and those Slater determinants whose energy is within 100 
meV of the lowest energy determinant are retained in the 
calculation. This rejection step reduces the total num- 
ber of determinants slightly and saves some computer 
time. The accuracy of the numerically calculated ground 
state addition energies is estimated to be 0.1 meV or bet- 
ter. The i?-dependent cut-off on n leads to small steps 
in the energy as a function of B which result from the 
sudden change in the number of basis states that occurs 
whenever the integer part of 6.9 — 0.285 changes. These 
artefacts are typically around 0.01 meV, about an order 
of magnitude smaller than experimental error. 

IV. DATA ANALYSIS 
A. Parameter fitting 

Standard techniques are used to fit the model parame- 
ters. The general idea is to find parameters that minimise 
the metric distance between experimental and theoreti- 
cal curves. This problem occurs in many areas of physics 
and the present work is based on a formalism developed 
for surface structure analysis [2l|, For each curve, 
the metric distance is taken to be 

^=]f E[/*(^'0-/e(i?.)f , (14) 

where ft and fe are theoretical and observed values of 
the addition energy or chemical potential at magnetic 
field Bi and Np is the number of experimental points 
in the curve. The metric distance used to fit multiple 
curves is the average of the metric distances for the in- 
dividual curves weighted with the number of points in 
each curve. With this choice, the procedure reduces to 
an unweighted least squares fit of the entire data set. It 
is possible to adjust the metric to give weight to specific 
features in the data, such as peaks and shoulders, but 
this is not done in the present work. The metric distance 
is minimised with the Levenberg-Marquardt algorithm, 
a standard numerical technique (23| . This fitting pro- 
cedure is generally accepted to give reliable parameter 
values provided the number of features in the data ex- 
ceeds the number of fitting parameters, a condition that 
is satisfied in the present case. 

B. Analysis of addition energy data 

The addition energy data was collected as part of an 
experiment to investigate ground state transitions in the 
high-magnetic field regime [H] ■ The data was taken with 
a fairly high magnetic field resolution (0.05 T) and con- 
sists of a limited number of VgN curves, = 2,3,4,5 
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FIG. 3: Experimental and theoretical addition energies. The 
field range used for parameter fitting is < B < 14 T. Dia- 
monds: experimental data, solid lines: theoretical results. 

only. The measurements were performed in a dilution 
refrigerator and the electron temperature is estimated to 
be below 100 mK. 

The ground state addition energy is the difference of 
two successive chemical potentials: Ean = En+i — 
2Em + En-1 — lJ.N+1 — IJ-N- It follows from Eq. ^ that 
Ean at fixed magnetic field is related to two successive 
gate voltages: 

Ean = e[a{VgN+i)VgN+i ~ a{VgN)VgN], (15) 

which is valid provided that the contact chemical poten- 
tial He is independent of Vg and therefore cancels when 
the aVg values are subtracted. Eq. (fT5|) can in principle 
be used to find Ean from experimental values of a and 
Vg . However this requires very accurate values of both a 
and Vg. Typical values of Ean are around 3 meV, while 
typical values of eaVg are around 60 meV. So both a and 
Vg need to be measured to an accuracy of 0.1 % or bet- 
ter to measure Ean to ±0.1 meV. Vg is known to this 
accuracy for the present device but a is not. Eq. (fT5|) is 
therefore approximated by 

Ean ^ eaN+i{VgN+i ~VgN), (16) 

where a.N+i = [a(V^Ar+i) + a(V^Ar)]/2. This suppresses 
the random errors which would result from use of Eq. (jlSp 
but is only valid when a(VgM) varies sufficiently slowly 
with N . This is assumed to be the case. The validity of 




IVlagnetic field (T) 

FIG. 4: Experimental and theoretical addition energies. As 
Fig. [3] except that the field range used for parameter fitting 
is < B < 10 T. 



this assumption cannot be tested without more data on 
a however the results based on this assumption are found 
to be consistent with results of chemical potential analy- 
sis, section PVCI Another factor that limits the present 
analysis is that a depends on the magnetic field as well as 
Vg but is only known for a sample of field points at inter- 
vals of around 1 — 2 T. The a values at intermediate field 
points have to be obtained by linear interpolation and 
this introduces a further systematic error which is difh- 
cult to quantify. In principle, both this systematic error 
and the one resulting from use of a can be eliminated by 
performing detailed measurements of a. 

Figure [3] shows the comparison between experimental 
and theoretical addition energies for A'' = 2, 3, 4. All 
three curves were used simultaneously to fit the param- 
eters, as described in section IIV Al It is clear that the 
absolute values of the experimental and theoretical data 
agree very well. The raw experimental data is shown in 
the figure together with raw theoretical data. There is no 
shifting or scaling of the theoretical curves. The agree- 
ment is generally good but there are some discrepancies, 
particularly around 4-6 T and above 10 T. To investigate 
the cause of these discrepancies the data was fitted again 
in the restricted field range < S < 10 T. The results 
are shown in Fig. [4] and it is clear that the fit is better, 
particularly around 4-6 T where the double peak struc- 
ture for A^ = 3 and A" = 4 is reproduced very well. The 
improvement in the fit can be quantified with the metric 
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distance which is equivalent to the RMS difference be- 
tween experiment and theory. The value corresponding 
to Fig.[3]is 0.23 meV while for Fig. Hit is 0.15 meV. 

The deterioration in the quality of the fit when the full 
data range is used could be caused by factors not included 
in the model or uncertainties in the experimental data. 
The most likely causes are impurity effects, screening and 
insufficient knowledge of a. 

Both impurity effects and screening could depend on 
magnetic field. Impurity effects are not included in the 
present model but are likely to be more significant in 
the very strong field regime. The mean distance between 
the impurities in the heavily doped contacts is around 
10 nm. This gives a fluctuating contribution to the dot 
potential on a similar length scale. But at zero magnetic 
field the typical diameter of the dot state is around 20 
- 30 nm. Hence the effect of the impurities will tend to 
average out. Indeed the experimental evidence is that 
the present device has a high degree of circular symme- 
try [5|. However as the magnetic field increases the dot 
wave function shrinks and eventually becomes smaller 
than the length scale of the potential fluctuations. Im- 
purity effects would be large in this regime but the mag- 
netic fields used in the present work are probably too low 
for this to happen. Magnetic field dependent screening 
is another possible cause of discrepancies but if the mag- 
netic field dependence was significant, it would probably 
occur throughout the field range, while the discrepancies 
are mainly in the high field range. 

The remaining possible cause of the discrepancies is 
that the value of a is not known to sufficient accuracy 
in the high-field regime. In this case the parameters de- 
termined by fitting the data up to 10 T are expected to 
provide an accurate description of the dot in the field 
range up to 14 T but the problem lies in the numerical 
value of the experimental addition energy. There is some 
evidence that this is the case: excited state features pre- 
dicted by the present model are found to correspond to 
features in experimental excitation spectra in the high- 
field regime Q . The value of a is not needed to identify 
these features and this suggests that the discrepancies 
in the high field addition energy are most likely to be 
related to insufficient knowledge of a. 

The best fit parameter values corresponding to Fig. [3] 
are a = 4.8 ± 0.1 meV, b = 0.02 ± 0.01 meV and = 
15.0 ± 3 nm. Here the statistical errors correspond to 
the parameter range compatible with the fiuctuations in 
the data shown in Fig. 3] but exclude the unquantifiable 
uncertainty in a. These values are taken to give the best 
model of the present device. The value of b suggests that 
the confinement energy increases slowly with gate voltage 
but the opposite would be expected from electrostatics 
of the device. However the value of b is comparable to 
its error so a model with constant confinement or slowly 
decreasing confinement would probably give an equally 
good fit to the data. 
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FIG. 5: Experimental and theoretical chemical potentials. Di- 
amonds: experimental data, solid lines: theoretical results. 

C. Analysis of chemical potential data 

In principle, Eq. ^ can be used to obtain iin from the 
Vg data but the appearance of the contact chemical po- 
tential in this equation presents an obstacle. The abso- 
lute value of fic is not known. In addition, /i^ is expected 
to vary with magnetic field and its field dependence is not 
known. Part of the difficulty can be eliminated by taking 
the chemical potential relative to the chemical potential 
for the first electron at zero magnetic field. The contact 
chemical potential is written as /^c(0) -I- Afic{B), where 
A/ic vanishes when B = Q. Eq. ^ becomes 

e[a{VgN, B)VgN{B) ~ a{VguO)Vgim = 

^iN{B) - ~ ^^i,{B). (17) 

This enables experimental and theoretical values of 
fiN{B) — /ii(0) to be compared provided that a model for 
A^c{B) is available. A model investigated in the present 
work is ApLc{B) = ^T'^ + (fiwc/2)2 - F, where F is a 
fitting parameter. An advantage of analysing chemical 
potential data is that Eq. (fT7|) is less sensitive to errors 
in a than Eq. (fT5|) however the disadvantage is that the 
results depend on the model chosen for A^c{B). 

Chemical potential data analysis requires a data set 
that contains Vgi so that fiN{B) — /xi(0) can be found. 
The addition energy data analysed in section IIVBI does 
not contain = 1 data so a second data set is used 
for the chemical potential data analysis. This data was 
collected under similar conditions to the addition energy 
data but over a wider range of N, a smaller magnetic field 
range and a lower field resolution (0.1 T). In addition, the 
device was subjected to one thermal cycle between col- 
lecting the two data sets. The chemical potential data 
was collected first then the device was taken up to room 
temperature and cooled again to collect the addition en- 
ergy data. 

The comparison between experimental and theoretical 
chemical potentials is shown in Fig. \5\ As with the ad- 
dition energy, a very good fit is obtained without any 
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scaling or shifting of the theoretical data. The RMS dif- 
ference between the experimental and theoretical curves 
is about 0.12 meV. The best fit parameter values are 
a = 4.8 ± 0.1 meV, b = -0.06 ± 0.02 meV, = 16.0 ± 2 
nm and F = 7.4±0.3 meV. Except for the sign of b, these 
values are consistent with the values obtained by fitting 
the addition energy. As with the addition energy the ab- 
solute magnitude of b is small and the statistical error in 
6 is relatively large. This suggests that the variation of 
hujQ with N is not very significant for the small N range 
considered here. 

The best model of the device is believed to be the 
model that results from fitting the addition energy be- 
cause an extra parameter is needed to fit the chemical 
potential data and it is preferable to keep the number of 
fitting parameters to a minimum. The addition energy 
model was therefore used in ref. Q to determine ground 
state quantum numbers from transport data. However 
ref. [5| contains excited state data and comparison of 
features in experimental and theoretical excitation data 
requires comparison of chemical potential data. This ap- 
pears to be an insurmountable problem because the addi- 
tion energy data does not contain Vgi{B) so it is impos- 
sible to fit r for the experimental conditions of ref. Q . 
However, the exact value of the chemical potential is not 
needed for the analysis of ref. [3| and similar situations 
because the only information that needs to be extracted 
from the experimental and theoretical data is the mag- 
netic fields at which features occur. As an aid to identi- 
fication of these features it is very useful to compute the 
chemical potential as hn{B) — /xi(0) — Aij,c{B) because 
this ensures that the experimental and theoretical data 
have roughly the same slope. This can be achieved by 
using an approximate value of F that is found by visual 
comparison of the data, 8.2 meV in the case of ref. @. 
The precise value of F is not important in this case as F 
does not affect the position of features in the transport 
data. 



V. APPLICATIONS OF THE MODEL 
A. Low field data 

The model parameters are found by fitting to the 
ground state addition energies so comparison of theoret- 
ical and experimental excitation energies is a good test 
of the predictive power of the model. This comparison is 
performed here for the low magnetic field regime, B < 5 
T, to which the model has not been applied before. A 
discussion of the high field regime is available in the lit- 
erature fe], see also section IVBl 

The experimental excitation energy data Fig. [5] (left 
panel) consists of the derivative of the source-drain cur- 
rent, dIsd/dVg, plotted as an intensity image in the B, Vg 
plane. The electron number ranges from 1 to 6. The data 
was collected in a similar way to the addition energy data 
but during a different thermal cycle. The source-drain 



voltage Vsd = 2 mV so the maximum excitation energy 
that can be probed is 2 meV. For each electron number 
there is a stripe that corresponds to this 2 meV window. 
The width of each stripe depends on B and Vg because of 
the B and Vg dependence of a. Within each strip there 
are dark lines that correspond to transport through ex- 
cited dot states. The excitation energy is found by scaling 
the line position to the stripe width. This gives excita- 
tion energies, AiJ, to an accuracy of about ±0.1 meV, 
together with field values accurate to about ±0.02 T. In 
addition to the excitation lines, there is 'noise'. Some of 
this is caused by emitter states which result from density 
of states fiuctuations in the heavily doped contacts. To 
extract the B dependence of the excitation lines it is nec- 
essary to distinguish contact features from dot features. 
This can be done by rejecting features that do not de- 
pend on A'^ for example. However the excitation lines 
in Fig. [S] have gaps and are of variable intensity so the 
present quantitative comparison is based on features like 
level crossings and crossings of excitation lines with the 
upper boundaries of the excitation stripes. These fea- 
tures can be identified unambiguously by comparing the 
topology of the experimental and theoretical data. 

The theoretical excitation energy data is shown in the 
right frame of Fig. El The figure shows ^n{B) — Aii(O) — 
A/ic(i?) for electron numbers in the range 2 to 5. For 
each N the energies of the ground state and lowest four 
excited states are calculated and the value of ij,n{B) — 
/ii(0) — Afj,c{B) is found for the ground state and each 
of the excited states. The excited state lines that fit into 
the 2 meV excitation window are shown by the solid lines 
in the figure. The continuous solid line at the bottom of 
each stripe is the ground state chemical potential and 
this coincides with the lower boundary of each excitation 
window. The dashed line at the top of each stripe shows 
the upper boundary, 2 meV above the ground state line. 

Qualitatively, there is very good agreement between 
the theory and experiment. However some of the exper- 
imental excited state lines are not continuous and there 
is some noise. The singlet-triplet transition of the 2 elec- 
tron system is clearly visible (a) and so is the excited 
state triplet line (b-a). For 3 electrons there are ground 
state transitions in the 4-5 T interval in both the the- 
ory and experiment and an excited state line is visible 
in the experimental data. This probably corresponds to 
line (c) in the theory. In the case of 4 electrons there 
are ground state transitions (d) and (e) in both theory 
and experiment. In addition, an excited state line which 
has a fiat maximum is present in the experimental data. 
This corresponds to the crossing (h) in the theory but 
small symmetry breaking effects such as disorder change 
the crossing into an anticrossing and lead to a maximum 
in the experimental data. The second excited state line 
(f-g) is also clearly visible in both the theory and experi- 
ment. For 5 electrons the ground state transition (j) has 
the distinctive form of a maximum and is clearly present 
in both theory and experiment, and some short excited 
state lines emerge from this crossing. Further, for both 
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FIG. 6: Experimental and theoretical excitation spectra. Left frame: experimental excitation spectra. The grey scale image 
shows dIad/dVg in the B, Vg plane for A'^ = 1-6. Right frame: theoretical excitation spectra for = 2-5. The dashed lines show 
the upper bound of the excitation stripe, AE — 2 meV. Roman font labels indicate the features listed in Table Hill 



TABLE III: Comparison of features in experimental and theoretical spectra. Features (a), (d), (e) and (j) are ground state 
transitions. For these features {L,S) — > {L',S') indicates low field values of {L,S) followed by high field values. Features (b), 
(f), (g) and (h) are excitation energies at the points shown in Fig. |S] and defined in the text. Feature (c') is the excitation 
energy on line (c) at _B = 3.0 T. For these features (L, 5") (L', S') indicates excited state {L, S) values followed by the ground 
state values. 



Feature N AEe^p (mcV) AExhe (meV) Bsxp (T) Brhe (T) (L, S) values 

(a) 2 0.0 0.0 3.8 4.2 (0,0)^(1,1) 

(b) 2 2.0 2.0 0.75 1.15 (1,1)^(0,0) 
(c') 3 1.4 1.22 3.0 3.0 (3, |) ^ (1, i) 

(d) 4 0.0 0.0 0.33 0.3 (0,1)^(2,0) 

(e) 4 0.0 0.0 4.06 4.0 (2,0) -^(3,1) 

(f) 4 1.0 1.24 0.0 0.0 (0,0)=>(0, 1) 

(g) 4 2.0 2.0 0.86 0.70 (0,0)=>(2,0) 

(h) 4 1.5 1.71 1.57 1.3 (0,1)^(3,1) 
(j) 5 0.0 0.0 1.39 1.25 (1, i) (4, i) 



4 and 5 electrons the experimental ground state line is 
broadened in the regions between 4 and 5 T. This cor- 
responds to the region where the theory predicts that 
ground and excited state levels cluster together and exci- 
tation gaps shrink. Although the agreement is generally 
good, there are some discrepancies. For TV = 1 the the- 
ory predicts no excitations within the 2 meV window at 
fields below 5 T but there is some structure in the exper- 
imental data. This is thought to be caused by emitter 
states. Similarly, for A'^ = 5 no excitations are predicted 
in the 2 meV window at zero magnetic field but there is 
some structure in the data. 

The quantitative experiment-theory comparison for 
features (a)-(j) is detailed in Table IIIII For each fea- 
ture, experimental and theoretical values of AE and B 
are given but AE = in the case of ground and excited 
state transitions. The theoretical magnetic fields are ac- 
curate to about 0.05 T, the magnetic field step used for 
the calculations. The theoretical quantum numbers are 



also listed. In the case of excitations (— >), the quantum 
numbers of the excited state are shown followed by those 
of the ground state. In the case of ground state transi- 
tions (=>), the quantum numbers on the low field side of 
the transition are followed by those on the high field side. 

In the case of the singlet-triplet transition of the 2- 
electron system (a) , the predicted transition field is about 
0.4 T lower than in the experiment and the point where 
the triplet excitation line crosses the 2 meV excitation 
boundary (b) is also about 0.4 T lower. This discrep- 
ancy is thought to be a consequence of the thermal cy- 
cling since the positions of the singlet-triplet transition 
agree well in the addition energy data, Figd) The excited 
state line of the 3-electron system (c) is difficult to iden- 
tify because it is not clear whether it extends as far as the 
excitation boundary in the experimental data. However 
there is some weak excited state structure in the data 
at around 3.5 T and this suggests that the experimen- 
tal line is the second excited state. The experimental 
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excitation energy is compared with the theoretical sec- 
ond excited state energy at the arbitrarily chosen field 
of 3 T in Table Hn] (feature (c')). The two values dif- 
fer by only 0.18 meV and this supports the idea that 
the experimental line is the second excited state. In ad- 
dition, the extrapolation of the experimental line crosses 
the upper boundary of the excitation stripe at about 2.24 
T compared with 2.4 T in the theory and line (c) crosses 
the theoretical ground state line at 4.6 T close to an ex- 
perimental ground state transition at 4.58 T. This gives 
further support to the interpretation of the experimental 
excited state line as the second excited state line corre- 
sponding to theoretical line (c). 

The 4-electron system is rich in features that can be 
identified unambiguously and compared quantitatively 
with theory. Ground state transitions at 0.33 and 4.06 
T agree very well with transitions at 0.3 T (d) and 4.0 
T (e) in the theory. The excitation energy of the second 
excited state (f ) at zero magnetic field agrees with theory 
to 0.24 meV while the position of the crossing of the sec- 
ond excited state line with the excitation stripe boundary 
(g) agrees with theory to 0.16T. The excitation energy at 
the excited state crossing (h) agrees to 0.21 meV and its 
position agrees to 0.27 T. It is more difficult to identify 
features in the 5-electron data as there is more noise. 
However the position of the ground state transition (j) 
agrees to 0.14 T. 



B. High field data 

Application of the present model to data in the high 
magnetic field regime beyond the MDD has been dis- 
cussed in ref. 5]. However the fitting procedures de- 
scribed here are not detailed in that work. Here it is 
emphasised that identical procedures were used in the 
present work and the work described in ref. Q and that 
identical model parameters were used in both cases. One 
of the main findings of ref. Q is the existence of inter- 
mediate spin states in the regime beyond the MDD. The 
MDD is spin polarised but as the magnetic field increases 
a transition to a partly polarised state occurs, followed by 
a spin polarised state, then another partly polarised state 
and so on. The quantum numbers found for the N = 5 
dot studied in ref. Q are consistent with general ideas 
about symmetry and electron molecular states 0, Q ■ In 
the specific case of the = 5 partly polarised states, 
calculations of the pair correlation function indicate that 
there is a superposition of 4- and 5-fold symmetry with 
a dominant 5-fold component. Further details are given 
inrefs. [Hill. 

VI. CONCLUSION 

An accurate model of a vertical pillar quantum dot has 
been developed which is able to reproduce addition en- 
ergy and chemical potential data for individual quantum 



dots. Experimental addition energies are reproduced to 
an accuracy of around 0.15 meV and this is accurate 
enough to enable ground state spin and orbital angu- 
lar momentum quantum numbers to be determined from 
transport data. Although the three-dimensional device 
structure is accounted for, the resulting physical model 
reduces to one in which electron motion is restricted to 
the two lateral dimensions, the confinement is parabolic 
and the interaction potential differs significantly from the 
bare Coulomb potential. 

The model only contains 3 adjustable parameters. Of 
these the parameter that describes the A^ dependence 
of the confinement energy does not appear to be signif- 
icant over the small N range considered here and it is 
likely that a two parameter model would give similar ac- 
curacy to the present one. For dots with larger A^ it 
would be possible to include non-parabolic confinement. 
This would increase the number of model parameters but 
with larger A^ there would be more data so the additional 
parameters could probably be determined reliably. 

Although the model parameters are determined by fit- 
ting ground state data, the model is able to give a good 
description of the low lying excited states. The model en- 
ergies agree with experiment to about 10-20% and the po- 
sitions of features in magnetic field dependent data agree 
to similar accuracy. One of the factors limiting the ac- 
curacy of the present analysis is uncertainty in the value 
of the electrostatic leverage factor a. Accurate measure- 
ments of a would enhance the scope and applicability of 
the present analysis. 

The general approach described here is not limited to 
vertical pillar dots. In principle, similar analysis proce- 
dures could be developed for any type of dot, provided a 
good physical model is available. 
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APPENDIX A: NUMERICAL CALCULATION OF 
GREEN'S FUNCTION 

The functions / and g and the Wronskian W (Eq. fTU]) 
are needed to find the Green's function. To find / and g 
the system is divided into thin slices, such that e and qo 
are constant in each slice. In the present case this simply 
means that each layer of the structure is sub-divided into 
thin slices but the same numerical procedure is applica- 



12 



ble to any form of e{z). The nth sHce occupies the region 
between z„^i and z„ and is characterised by a relative 
permittivity e„, q value = 9^ + qo{z) and a thickness 
tn = Zn — Zn-i- Within each slice / and g have the form 
A„exp(— (7„z) + i3„exp(g„z). The two terms in this ex- 
pression are analogous to evanescent waves in diffraction 
theory and it is convenient to use the language of diffrac- 
tion theory to describe them. Thus the reflected wave, 
R which decays in the positive z direction has the form 
R = Anexp(—qnz). Similarly, the transmitted wave, T 
which decays in the negative z direction has the form 
T = B„cxp(g„z). 

The amplitudes at successive slices are related via a 
transfer matrix. 



R71+1 

Tn+l 



mi m2 
m3 7714 



Rn 
T 



(Al) 



where i?„ and r„ are amplitudes just below the inter- 
face at z„, that is within the slice of relative permittivity 
e„. The boundary conditions on G are used to find the 
elements of the transfer matrix: g and e{z)dg/dz are con- 
tinuous and the same holds for /. Therefore 



mi = - 1 



^ 1- 



7712 



"^3 = 7: 1 



"^4 = 77 1 



^nqn 

Cn+lgn+l 

^nqn 
Cn+lgn+l 



cxp(-9„+ii„+i), 
exp(-gf„+ii„+i), 
exp((7„+it„+i), 
exp(g„+it„+i). 



(A2) 



It is well known that direct application of Eq. (|Aip is 
numerically unstable. Instead it is necessary to compute 
the ratio r„ = Rn/Tn which corresponds to the reflection 
coefficient at each interface [13, lisl . ITgl . The notation 
used here is similar to that of ref. '19'|. The reflection 
coefficient r„ satisfies 



Tn+l 



mir„ 



TO2 



TO3r„ -I- TO4 



(A3) 



where the are the elements of the transfer matrix 
(Eq. (|A2[) ) for going across the interface at z„ to the in- 
terface at Zn+i- Eq. (jA3[) can be used to step r„ provided 
that the denominator of the fraction is not small. If the 
denominator does become small an alternative relation 
can be used to step [l3| but this was not needed for 
the present work. 

The procedure for finding g is as follows. Deep in- 
side the bottom contact g decays exponentially so the 
only component present is T„. Hence r„ is found from 
Eq. (jA3[) starting from tq = 0. Once Tn is known gn is 
found from the relations 



Tn = (TO3r„+TO4) ^Tn + l, 
= {l+rn)Tn, 



(A4) 



which follow from Eq. (|A1|1 and the definitions of i?„ and 
Tn- Eq. (|A4p is applied with the initial condition Tn — 
1, where the topmost slice ends at zat. This procedure 
corresponds directly to the one used in diffraction theory. 

The procedure used to find / is slightly different. / 
is required to decay exponentially deep inside the top 
contact, that is, it only has a reflected component there. 
This means that the inverse reflection coefficient T„/i?„ 
vanishes deep inside the top contact. Hence it is conve- 
nient to re-expresss Eqs. (jA3p and (jA4|) in terms of the 
inverse transfer matrix defined by 



Rn 
Tn 



nil n2 



Rn+l 
Tn+l 



(A5) 



where i?„ and r„ are amplitudes just below the interface 
at Zn-i, again within the slice of relative permittivity e„. 
The boundary conditions on G lead to expressions for the 
inverse transfer matrix elements: 



ni = i ( 1 
n2 = i I 1 



1 



n4 = - 1 1 



eri+lgiiH 


-1 


^nqn 




Cn+lgnH 


-1 


^nqn 




Cri+lgriH 


-1 


^■nqn 




en+lgnH 


-1 





cxp(g„+it„+i), 
cxp(-g„+ii„+i), 
exp(g„+it„+i), 
exp(-(7„+ii„+i). (A6) 



The inverse reflection coefficient satisfies 



113 + n4f„|i 



(A7) 



and this relation is used to step f~^, starting from the 
initial condition rj;^^i = 0. Here the notation r is used 
to emphasise that r and r are computed with different 
initial conditions so is not the inverse of r. Once 
is known, /„ is found from the relations 



Rn+i = {ni+n2r^l^) ^i?„, 

= (i + r;;i)i?„. 



(A8) 



with the initial condition Rq — 1. 

The Wronskian W is independent of z and can be eval- 
uated at any convenient position. This requires /, g and 
their derivatives. The derivatives are found from 



c5f 

dz 



qni-An exp(-g„z) + Bn exp((j„z)) 
qn{T-R) 



(A9) 



together with a similar expression for dg/dz. This leads 
to 

W = -2qe^e. ^^"["''""^ J n9n, (AlO) 
(l+r„ )(l + r„) 
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which is vahd away from dielectric interfaces. 

The relations given here have been found to be nu- 
merically stable for the device considered in the present 
work. But in general it is possible for the calculation 
of the Wronskian to underflow, causing overflow in the 
calculation of the Green's function. However this is only 
likely to be a problem for very thick systems and where 
underflow could result from the exponential form of / 
and g. 



APPENDIX B: NUMERICAL CALCULATION OF 
GREEN'S FUNCTION INTEGRALS 

The Fourier components of the effective interaction are 
found from Eq. pT|) . The form factor is 

F{q) - 2ge^eo j x^{z)x^{z')G{q, z, z')dzdz' . (Bl) 

Except for a factor of —1/W, the integral in this equation 
has the form 

X^{z).f{z) j x^{z')g{z')dz'dz 

■OO J — OO 

+ I x\z)g{z) / x\z')f{z')dz'dz 

J — OO J z 

/OO pZ 
X'{z)f{z) / x'{z')g{z')dz'dz. (B2) 
-OO J —OO 

To evaluate this integral numerically it is convenient to 
define 



and the required integral is found from 



xHz)I{z)dz. 



(B5) 



The advantage of this approach is that only a few opera- 
tions are needed to update I{z). As a result the computer 
time scaling is linear with the number of integration steps 
instead of quadratic as would be the case if the integral 
in Eq. (|B2p was evaluated directly. 



The approach is simple to implement on a uniform 
grid. However a uniform grid cannot be used in the 
present case because the system contains interfaces where 
the relative permittivity changes abruptly and it is neces- 
sary to ensure that the grid points coincide with these in- 
terfaces. This is done by adjusting the step length within 
each region of constant relative permittivity. As a result 
the step length is constant within each region but the 
step length varies from region to region. To apply the 
approach to the resulting non-uniform grid the following 
generalisation of Simpson's rule is used: 

2Q+A2 



u{z)dz ^ au{zQ — Ai) + fcu(zo) -I- cu{zq + A2), 

(B6) 

where a, 6, c are chosen so that the integration rule is 
exact for u{z) — 1, z, and Ai and A2 are step lengths. 
The values of a.b.c are 



I{z) = f{z) I x\z')g(z')dz' 



Then I{z) is evaluated recursively from 



f{z + Az) 



f{z) 
\\z')g{z')dz' 



(B3) 



(B4) 



2A2 + A1A2 



2Ai 



^2 



6A1 

3(Ai + A2)AiA2 



A? 



6A1A2 
A1A2 - A? 



6A2 



(B7) 
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